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Abstract 

We study the behavior of nonlinear waves in a two-dimensional medium with density and 
stress relation that vary periodically in space. Efficient approximate Riemann solvers are devel- 
oped for the corresponding variable-coefficient first-order hyperbolic system. We present direct 
numerical simulations of this multiscale problem, focused on the propagation of a single localized 
perturbation in media with strongly varying impedance. For the conditions studied, we find little 
evidence of shock formation. Instead, solutions consist primarily of solitary waves. These solitary 
waves are observed to be stable over long times and to interact in a manner approximately like 
solitons. The system considered has no dispersive terms; these solitary waves arise due to the 
material heterogeneity, which leads to strong reflections and effective dispersion. 

1 Introduction 

Solutions of first order hyperbolic partial differential equations (PDEs) generically lead to forma- 
tion of shock singularities and subsequent entropy decay. In contrast, nonlinear wave equations 
with dispersive terms, such as the Korteweg-deVries equation (KdV), may exhibit solitary wave 
solutions |14| . These solutions are remarkable in the context of nonlinear PDEs since they prop- 
agate without changing shape and, in many cases, interact only through a phase shift. Thus the 
long-time solution behavior of nonlinear waves depends critically on the presence of dispersive 
regularizing terms. 

Santosa and Symes [13 showed that solutions of the linear wave equation in a medium with 
periodically varying coefficients exhibit dispersive behavior, even though the equation contains no 
dispersive terms. LeVeque and Yong [11] found that solutions of the ID p-system in a periodic 
layered medium form solitary waves (referred to as stegotons due to their discontinuous shape) 
similar to those arising in nonlinear dispersive wave equations like KdV. This is remarkable since 
the equations they considered are first order hyperbolic PDEs with no dispersive terms. Recently 
it has been observed that solutions of first order hyperbolic systems with periodic coefficients 
may be characterized by shocks or by solitary waves, depending on the medium and the initial 
conditions [5 . Further experiments in [4j [5] indicate that solitary waves may arise quite generally 
in the solution of nonlinear hyperbolic PDEs with periodically varying coefficients. Those results 
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have revealed new behaviors of nonlinear waves in ID heterogeneous materials, and it is natural 
to ask if such behaviors persist in higher dimensions. The present work, which builds on [12 , 
extends the aforementioned studies to a fully two-dimensional setting. 

We consider nonlinear wave propagation in a 2D periodic medium, modeled by variable- 
coefficient first order hyperbolic PDEs. Since ID solitary waves are often found to be unstable 
when extended to higher dimensions, a principal question is whether waves like stegotons are 
stable in 2D, and whether they form in materials that vary in multiple spatial directions. In order 
to investigate these questions in a general setting, we have chosen to study the nonlinear wave 
equation 

e "- v -(^ v<T(e ' x) ) =0 - 

This is perhaps the simplest multidimensional nonlinear wave model that is general in the sense of 
allowing for wave propagation in all directions. In ID, it is commonly referred to as the p- system 
due to its connection with Lagrangian gas dynamics. 

In order to introduce the object of our study, four representative examples of the behavior 



of cylindrical wavefronts are depicted in Figure 1(c) Each solution shown results from the same 



initial condition: a small, cylindrically symmetric Gaussian pulse centered at the origin, shown 



(close-up) in Figure 1(a) and given by: 



/ m r ( (x-0.25) 2 (?/ - 0.25) 2 \ 
(7(a5,0) = 5exp f ^-L- - Ky 1Q ; j . (1) 

Since the solution is symmetric under reflection about the x- and ?/-axes, a single quadrant is 



sufficient to characterize the full solution. The top left quadrant of Figure 1(c) shows the solution 
obtained in a linear, homogeneous medium: a smooth pulse expanding at the sound speed of the 
medium. The top right quadrant corresponds to a nonlinear, homogeneous medium; the leading 
edge of the pulse steepens into a shock wave. The bottom right quadrant corresponds to a linear, 



heterogeneous medium shown in Figure 1(b) that is composed of alternating square homogeneous 
regions in a checkerboard pattern. The front travels more slowly than in the homogeneous case, 
due to the effect of reflections. Finally, the bottom left quadrant shows the subject of the present 
study. The medium is both nonlinear and heterogeneous (with the same checkerboard structure) . 
In this case, a train of cylindrical pulses forms. As we will see, these pulses appear to behave as 
solitary waves and we refer to them as cylindrical stegotons. Note that they do not exhibit exact 
cylindrical symmetry, due to the lack of such symmetry in the medium. 

In the remainder of the paper, we investigate the conditions that lead to these cylindrical stego- 
tons, as well as their properties and behavior. The model equations and materials are introduced 
in Section |2J 

Accurate numerical modeling of nonlinear waves in multidimensional heterogeneous media is 
challenging. Care must be taken to properly handle material inhomogeneities (especially dis- 
continuities) and shocks. We use finite volume methods based on approximate Riemann solvers 
implemented in Clawpack [8] and SharpClaw [7], New approximate Riemann solvers for the het- 
erogeneous 2D p-system are discussed in Section [3] Because very large grids (with approximately 
7 billion unknowns) are required to resolve this multiscale problem over times and distances suffi- 
cient to observe solitary wave formation, we have used the parallel PyClaw framework [6] in order 
to run on 16,384 cores of the Shaheen system at KAUST. 

Motivated by known behaviors of more traditional solitary waves, and also by studies of ID 
stegotons, we investigate qualitatively the following questions: 

• What kinds of media and initial conditions give rise to solitary waves? 
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(a) Initial condition (close-up) 
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(b) Checkerboard medium 
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(c) Solution in four different media 



Figure 1: Wave behavior in different media. Top left: linear, homogeneous. Top right: nonlinear, 
homogeneous. Bottom right: linear, periodic. Bottom left: nonlinear, periodic. The quantity plotted 
is stress (a). 
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• How do cylindrical stegotons interact with each other? 

• What is the role of shock formation in 2D periodic media with strongly varying impedance? 
These questions are addressed through further numerical experiments in Section [4] 



2 The 2D spatially- varying ^-system 
2.1 Governing equations 

As mentioned in the introduction, we have chosen to study the p-system due to its relative 
simplicity and generality. This system can be viewed as a simplification of models governing more 
complex wave behavior, such as that of elastic waves. In a solid, an elastic wave is composed 
of longitudinal or P- waves and transversal or S- waves [9]. If we assume the stress is hydrostatic 
(i.e., there is no shear stress and the extensional stress components are equal), the propagation 
of elastic waves can be modeled by: 



€tt 
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1 



-Vcr(e,x; 



0, 



(2) 



where e represents the strain, p(x) is the spatially- varying material density, <j(e, x) is the stress 
and x = [x,y] T is the position vector. Similar to [11] , we consider the nonlinear constitutive 
relation 



cr(e, x) = exp(if(x)e) + 1, 



(3) 



where K(x) plays the role of bulk modulus. Equation ^ with the stress relation ^ admits 
shock formation. In order to determine entropy-satisfying weak solutions, let us write (J2| as a 
first-order hyperbolic system of conservation laws: 



Qt + f (q, x)^ + g(q, x) y = 0, 



(4a) 
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(4b) 



Here u and v are the x- and ^/-components of velocity, q is the vector of conserved quantities, 
and f , g are the components of the flux in the x- and ^/-directions, respectively. This form arises 
naturally through consideration of kinematics and Newton's second law, and leads to the correct 
jump conditions across shocks. 



2.2 Periodic Media 

We consider two types of periodic variation in the density p(x) and the bulk modulus K(x). The 
period of the medium is always taken to be unity. The first medium, shown in figure 2(a) consists 
of a checkerboard pattern: 
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(a) Checkerboard medium (equation j5j) 
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(b) Sinusoidal type medium (equation 
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Figure 2: The domains discussed in Section [272] Sections of size 10 x 10 are shown here, but much 
larger domains are used in the numerical simulations. 



The second medium, shown in figure 2(b) is similar but smoothly (sinusoidally) varying: 
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3 Numerical Discretization 

In this section we describe the numerical methods used to solve the 2D spatially varying p- 
system Q. For all computations in this work we use the algorithms implemented in Clawpack 
and SharpClaw. Both are finite volume methods for solving hyperbolic PDEs based on solving 
Riemann problems. For the simulations presented here, both Clawpack and SharpClaw are driven 
through PyClaw [6] which is a lightweight Python framework that calls the low-level Fortran 
routines of Clawpack and SharpClaw, and also interfaces with PETSc [1 to provide parallelism. 

After briefly introducing the discretization schemes, we focus on their key ingredient, the 
approximate solution of the Riemann problem for the spatially- varying p-system Q . The Riemann 
solvers developed here follow the ideas of [10]. In particular, they are based on use of an all-shock 
solution and /-wave decomposition [2]. 



3.1 Second-order Clawpack discretization 

Clawpack is based on Lax-Wendroff discretization combined with TVD limiters, and is second 
order accurate in space and time [8]. The multidimensional Clawpack algorithms used here 
require the propagation of waves in both the normal and transverse directions at each cell edge. 
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The Clawpack discretization takes the form 

At 



_ At 
Ay 

At f~ n ~ n \ At (~ n ~ n \ 

where ^4 ± AQ n i . and B ± AQ n i are first-order fluctuations computed by the normal Riemann 



solvers described in Section 



The quantities F n ,i . and G n ._i are second-order corrections 



3.3 

that include the fluctuations computed by the transverse Riemann solvers described in Section 



3.4 as well as high-resolution approximations to c\x X ,c\y y . For details, see [9]. 



3.2 High-order WENO discretization using SharpClaw 

SharpClaw is based on the method of lines approach and uses fifth-order weighted essentially 
non-oscillatory (WENO) reconstruction in space with fourth-order strong stability preserving 
(SSP) Runge-Kutta integration in time [7]. SharpClaw requires propagation of waves only in the 
direction normal to each edge. 

First, a WENO reconstruction of the solution is computed from the cell averages Qi to give 
high order accurate point values q L _i and q K _i just to the left and right (respectively) of each cell 

% 2 % 2 

interfaces x i _i. A Riemann solution is computed at each interface based on the reconstructed 
values there. The resulting fluctuations are used to update the adjacent cell averages. An addi- 
tional term appears that is proportional to f x t+ ^ Aq x dx. For conservative systems like Q, this 

l ~ 2 

term can be conveniently computed in terms of a fictitious internal Riemann problem in each cell 
[7]. The semi-discrete scheme takes the form 



n+1 



dt Ax V ^+2>^ ^-2>J ^ LJ J Ay 

SharpClaw employs the same wave propagation Riemann solvers and user interface as Claw- 
pack. In multi-dimensions SharpClaw requires propagation of waves only in the normal direction 
to each edge. 



3.3 Normal Riemann solvers 

We assume that the density and bulk modulus are constant within each grid cell: p = Pij,K = K i3 - . 
In the case of smoothly varying media, this is an approximation. Then the system of conservation 
laws p]) can be written in the following quasilinear form within cell 



q t + Aij<i x + Bijtiy 
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Note that a e ^j denotes the derivative of <j(e, Xi, yj) with respect to e. For concreteness, we consider 
a Riemann problem in the x-direction, at interface x i _ 1 ■ which consists of the hyperbolic system 
Q with coefficients 



and initial condition 



Pi- 1,3 
Pi, 3 



Vx < x { 
Vx > x- 



q(z,y) 



Qi- 

Qi,3 



1,3 



K(x,y) 



Vx < x-. 
Vx > x- 



1,3 



Vx < x- 
Vx > x- 



(9) 



(10) 



In practice, Qij may represent a cell-average (in Clawpack) or a reconstructed value at the cell 
interface (in SharpClaw). The eigenvectors of A^ (Q^) are 







with corresponding eigenvalues {— dj, 0, Here Zij 

is the sound speed. 



y/Pij(7e,ij is the impedance and q. 



In the linear case, each wave in the Riemann solution is a discontinuity proportional to the 
corresponding eigenvector in the material carrying the wave. Thus the solution can be found by 
decomposing the difference AQ^_ i = — Qi-i^ in terms of the following three eigenvectors: 



-1,3 



1 





-1' 





(11) 



In the nonlinear case, the Riemann solution also consists of three waves, one of which is a 
stationary shear wave with zero velocity. The other two waves may be rarefaction waves or shock 
waves, but the solution cannot include transonic rarefactions, since the p- system is derived in a 
Lagrangian frame. For reasons of computational efficiency, we will use an approximate all-shock 
solver. Shock waves in the Riemann solution correspond to traveling discontinuities that are pro- 
portional to (left-going) or rf ;J - (right-going). Equations for an exact all-shock Riemann 
solution can be derived in a straightforward manner by considering the Rankine-Hugoniot condi- 
tions. However, these lead to a coupled nonlinear system that is expensive to solve numerically. 
Instead, we again follow [10] and approximate the solution further by replacing the exact wave 
speeds with the sound speeds in each cell. The waves themselves are found following the /-wave 



approach [2], by decomposing the jump in the normal flux in terms of the eigenvectors ( 11 ) 



-1,3 



2<1 
3 



where = f (Qij). The second wave is not needed in the numerical solution, since it has velocity 
zero and thus does not affect the solution value in either cell. The other two waves are accumulated 
into left- and right-going fluctuations, which are quantities that consider the net effect of all left- 
and right-going waves respectively: 



(12) 
(13) 
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The normal Riemann solver for the y-direction uses the same approach. To solve the Riemann 



problem at in place of (11) we use the eigenvectors 



1 1 

r ^--i = r *,;-i = 
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The normal flux difference is decomposed in terms of these eigenvectors: 
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(14) 



where Gij = g(Qij)- Finally, the waves are accumulated into up- and down-going fluctuations: 

B"AQ.._ 1=^,1 (15) 



B + AQ r<J 



(16) 



As observed in [10] for the ID p-system, this approach leads to efficient approximate Riemann 
solvers that are conservative and provide good accuracy at least for weakly nonlinear problems. 

3.4 Transverse Riemann solvers 

The multidimensional Clawpack algorithm makes use of a transverse Riemann solver that de- 
composes the horizontally traveling waves into up- and down-going corrections and the vertical 
traveling waves into right- and left-going corrections. These second-order corrections capture the 
effect of corner transport. 

Transverse corrections to the vertical fluctuations are computed by decomposing the horizontal 



fluctuations into up- and down-going parts using the eigenvectors (14) 



A ± AQ i _ 1 
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(17) 



The speed of the waves in the vertical direction determines the amount of horizontal fluctuation 
that must be added to the vertical one. These speeds are given by the eigenvalues of B ij -_i, 
which are s 1 . 1 = —c- ■ 1 , s 2 . 1 = and s 3 . 1 = c • • 1 . Finally, the corrections to the 
vertical fluctuations are: 



B~A ± AQ i 1 
B + A ± AQ i _i 
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The transverse corrections to the horizontal fluctuations are obtained in an analogous way. 
Those corrections are: 

A~B ± AQ 1 = 4_i ,7-_i 

hJ 2 2 '"^ 2 2 

3 



^ + B ± AQ i , 1 = 7 ?_i r,_ . 
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3.5 Accuracy tests 



For the nonlinear, variable- coefficient system studied in the present work, exact solutions are not 
available and error and convergence estimates are difficult to obtain because the degree of regu- 
larity of solutions is not known. In this section some measure of the accuracy of the numerical 
solutions is obtained by conducting self- convergence tests on small problems with the same qual- 
itative features as the problems of interest. The principal purpose of these tests is to determine 
the grid resolution necessary to ensure small relative errors and qualitatively accurate results. 
Detailed numerical analysis of the schemes' convergence behavior for these problems is beyond 
the scope of this work. In all computations here and in the rest of this work, a uniform cartesian 
grid is used with Ax = Ay = h. 



3.5.1 Sinusoidal medium 

Considering the nonlinear constitutive relation Q and the material properties pB = Kb = 10, we 
perform a self-convergence study in the sinusoidal medium |6| using the Clawpack discretization. 
The initial condition is given by The computational domain is restricted to the positive 

quadrant by imposing reflecting boundary conditions at the left and bottom. The stress at t — 3 
on a grid with h* — 1/480 is used as reference solution. Table [I] shows self-convergence rates as 
well as relative errors, computed by: 

\\a-a*\\2 , v 

~h' Mb ' ( } 

where a is the computed solution and a* is the solution on the fine grid. The observed convergence 
rate is roughly second order. 
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L 2 Error 


Rate 
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1.737xl0- 3 
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8.943xl0~ 4 


1.638 


160 


5.306xl0- 4 


1.814 


240 


2.235xl0~ 4 


2.132 



Table 1: Self-convergence test for the sinusoidal medium using the Clawpack discretization. 



3.5.2 Checkerboard medium 



Next we consider self- convergence of the solution of Q in the checkerboard medium |5| using 
the Clawpack discretization and SharpClaw. We use the nonlinear constitutive relation given 
by |3| and material parameters pB = Kb = 5. The initial condition is given by |l]). The 
stress is computed at t = 3 and the solution on a grid with h = 1/480 is taken as the reference 
solution for the self-convergence study. Tables 1(a) and |l(b)| show the convergence rates and 



relative errors following (18) using Clawpack and SharpClaw respectively. The achieved order 
of convergence is between 1 and 2 for both discretizations, which seems reasonable given the 
discontinuous coefficients. The SharpClaw solution is noticeably more accurate, even though the 
two methods exhibit similar convergence rates. 
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(a) Clawpack 






(b) SharpClaw 
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Table 2: Self-convergence test for the checkerboard medium using (a) Clawpack and (b) SharpClaw. 



4 Computational results 
4.1 Formation of solitary waves 

We now investigate in detail the solitary wave trains mentioned in Section [I] These waves are 
found to arise in both the checkerboard and the sinusoidal domain. We take a symmetric Gaussian 
hump |l]) as initial stress and solve the p-system |2| with the nonlinear constitutive relation |3| . 
Reflecting boundary conditions are used at the left and bottom boundaries in order to take 
advantage of symmetry and compute only in the first quadrant. The mesh is uniform with 
h = . Based on the accuracy studies above and additional tests, we expect this to be sufficient 
computational resolution to obtain qualitatively accurate solutions. 

For the checkerboard domain, we take pB = Kb = 5 and compute the solution using SharpClaw 
with h = ^0 . For the sinusoidal domain, we take pb — Kb — 10 and compute the solution using 
Clawpack with h — Figure [3] shows the stress at t = 90 and t = 200 for both materials, 

including slices along the lines y — x and y = 0. The persistence of the solitary waves after long 
times strongly suggests that they are stable, attracting solutions. 

Note that the solitary wave front is noticeably non-circular and that the cross-sectional shape 
of the solitary wave pulse varies with respect to the angle of propagation. 



4.2 Entropy evolution and shock formation 

As discussed in the introduction, solutions of nonlinear hyperbolic PDEs generically develop shock 
discontinuities, leading to irreversibility and entropy decay. It has been observed that spatially 
varying materials can inhibit the formation of shocks [3] [5]. Meanwhile, dispersive nonlinear wave 
equations commonly exhibit solutions that are regular for all times. Here we are solving a PDE 
without dispersive terms; nevertheless, the spatially varying coefficients create reflections, which 
yield effective dispersion. 

No visible shocks appear in the solutions shown in Figure|3] On the other hand, figure [4] shows 
that if the impedance ratio is small enough, the solution develops shocks. This is in qualitative 
agreement with the theory for ID systems developed in [5]. We now study shock formation in the 
2D stegotons by considering the entropy evolution and the reversibility of the solution. 

An entropy function for a hyperbolic PDE is a function that is conserved while the solution is 
smooth, but decreases in time when shocks form [9]. Thus, entropy decay is a signature of shock 
formation. A suitable entropy function for the 2D p-system |2]) is the total energy: 

r](u,e) = J (^p{x)u 2 + J cr{s,x)ds S j dx. (19) 
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(a) Checkerboard medium, t = 90, h = 1/240 




(b) Checkerboard medium, t = 200, h = 1/240 
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(c) Sinusoidal medium, t = 90, h = 1/240 




(d) Sinusoidal medium, t = 200, h = 1/120 



Figure 3: Stress at t = 90 and t = 200 for solitary wave trains in the sinusoidal and checkerboard 
media. The slice plots (on the right) show results along the lines y = x (solid black line) and y = 
(dashed red line). 
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Figure 4: (a) Stress at t = 90 in a checkerboard medium with small impedance ratio: pB = Kb = 1.5. 
(b) Slice plots showing results along the lines y = x (solid black line) and y = (dashed red line). 
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Figure 5: (a) Stress at t = 80 in a homogeneous nonlinear medium, (b) Slice plots showing results 
along the lines y — x (solid black line) and y = (dashed red line). 



We now compare the solution in a homogeneous domain with p = K — 1 (shown in Figure [5| to 
that obtained in the sinusoidal domain with pb = Kb = 10 (shown in Figure 3(c) ). The difference 
in the behavior is clear: the heterogeneity introduces reflections that effectively yield dispersion, 
breaking the initial profile into solitary waves. For the homogeneous domain, no reflections are 
present; consequently, there is no dispersion to prevent shock formation. In Figure [6] we see 
the entropy evolution for both cases. For the simulation with homogeneous domain, the entropy 
starts to decrease as soon as the shock starts to happen. On the other hand, the entropy using 
the sinusoidal type medium remains almost constant. In both cases we use a resolution given by 
fi — J— 

11 240 • 

Indeed, in both cases, there is loss of entropy due to numerical dissipation, but this converges 
to zero, as the grid is refined. To better determine if there is shock formation in the sinusoidal 
case, we present, in figure 7(a)| the normalized entropy evolution considering different resolutions 



given by ^ = 80, 120, 160, 240. We see the entropy lost is indeed small. However, it is evident that 
the rate of entropy lost changes before t = 10 making the entropy lost not to converge to zero as 
the grid is refined. This suggests the existence of a shock. To corroborate this, we focus up to 



12 




10 20 30 40 50 60 70 80 
t 



Figure 6: Entropy evolution considering a homogeneous nonlinear medium (solid line) and a het- 
erogeneous nonlinear medium (dashed line). The entropy has been normalized relative to the initial 
entropy. 



t = 10, study the entropy for even finer grids and see its convergence. In figure 7(b) we show the 
entropy evolution up to t = 10 for different resolutions given by \ = 80,120,160,240,480,960. 
We see the entropy at t — 10 converges to a value different than the initial entropy; i.e., there is 
a loss of entropy, which corroborates there existence of a shock. 

4.3 Interaction of cylindrical stegotons 



In order to study the interaction of these waves, we start with the solution depicted in Figure 3(b) 
We attempt to isolate the leading pulse in that solution by setting the solution to zero outside of 
a narrow band. Let q denote the solution of Q and q denote the corresponding isolated solution. 
Then q is given by: 

<\{xi, yj ) = l (20) 
[0 if Xi < X(j). 

Here X(j) is chosen so as to separate the leading stegoton from everything to the left as well as 
possible. Specifically, for each j, X(j) is the right-most local minimum of q: 

X(j) = max{x; : Sij < Si+ij;Sij < Si-ij}, (21) 

i 

where S is the discretization of the stress a. This isolation is imperfect since it seems clear that 
the tails of the solitary waves still overlap. Obtaining completely separated solitary waves would 
require even greater computational resources or a different modeling approach. By extracting just 
the leading pulse in this manner at two different times (t — 180 and t — 190), we obtain a pair of 
nearly isolated solitary waves. The interaction simulation is initialized using the sum of these two 
solutions as initial condition, but with the velocity negated in the t — 190 solution so that it will 



propagate inward. Figure 8(a) shows surface plots of the solution at t — 0, 3, 6, 9. Corresponding 



slices at y = x are shown in Figure 8(b) (solid line). For comparison, we also simulate just the 
outer pulse (without the inner pulse) over the same time (dotted line). 

The results are typical of solitary waves; the two pulses retain their identity after the interac- 
tion. Furthermore, the position of the outward-going pulse after the interaction is essentially the 
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(a) Entropy evolution up to t = 80 
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(b) Entropy evolution up to t = 10 



lire 7: Entropy evolution up to (a) t = 80 and (b) t = 10 considering (a) ^ = 80, 120, 160, 240 and 
j- = 80, 120, 160, 240, 480, 960. The entropy has been normalized relative to the initial entropy. 



same as in the simulation with no interaction. The apparent lack of a noticeable phase shift seems 
to be a result of the very short interaction time for this head-on collision. The same effect can be 
seen in one-dimensional stegoton interactions, as shown in Figure [9] Two stegotons traveling in 
the same direction exhibit a phase shift after interaction, but no phase shift is discernible after a 
head-on interaction. In order to simulate the interaction of two cylindrical stegotons traveling in 
the same direction (and the associated phase shift), much greater computational resources would 
be required. 

Some small oscillations are visible in the final solution in Figure [8] Since similar oscillations 
are seen in the solution obtained with just the outward- going pulse, it seems that these may be 
attributed to the fact that the tails of the solitary waves were not accurately captured in the 
initial condition. However, it might also be that some small oscillations are radiated during the 
interaction. 



5 Conclusion 

We have performed direct numerical simulations of the multiscale behavior of nonlinear waves 
in non-dispersive periodic materials. The largest simulations involve 6.9 x 10 9 unknowns and 
were performed on 16K cores of a BlueGene/P system. Our computational results indicate that 
solitary waves may arise in solutions of two-dimensional first-order nonlinear hyperbolic PDEs 
with spatially periodic coefficients, including both smooth and piecewise-constant media. As in 
ID, these waves apparently result from the combination of effective (material) dispersion and 
nonlinear steepening. The effective dispersion is a macroscopic effect caused by reflections in 
variable-impedance media. In media with strongly varying impedance, shock formation is strongly 
suppressed relative to that occurring in homogeneous media, as demonstrated by near-conservation 
of entropy. The cylindrical solitary waves that we have studied appear to behave approximately 
like solitons in their interactions. 

Studying the properties of these waves in detail is difficult due to the multiscale nature of the 
problem. Future work may include simulations using massively parallel adaptive mesh refinement 
or more computational resources in order to more fully isolate the solitary waves and study 
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(a) Spatial distribution of stress. 
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(b) Slices at y = x. 

Figure 8: 2D collision of cylindrical stegotons for t = 0, 3, 6, 9. The stegotons are initially moving 
toward each other. 
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(a) ID collision with same direction. 
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(b) ID collision with opposite direction. 
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Figure 9: Interaction of ID stegotons. The dashed line shows the stegoton originally at the left 
propagating on its own. 



longer interactions. A complementary tool in understanding these waves is the application of 
homogenization theory to multidimensional nonlinear hyperbolic systems. We have focused on 
localized perturbations in media with uniform sound speed and strongly varying impedance. Many 
other geometries, for both the medium and the perturbation, could also be considered. All of these 
topics are the subject of ongoing work. 



References 

[1] Satish Balay, Jed Brown, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G. 
Knepley, Lois Curfman Mclnnes, Barry F. Smith, and Hong Zhang. PETSc Web page, 2012. 
http:/ /www. mcs.anl.gov/petsc. 

[2] D.S. Bale, R.J. LeVeque, S. Mitran, and J. A. Rossmanith. A wave propagation method for 
conservation laws and balance laws with spatially varying flux functions. SI AM Journal on 
Scientific Computing, 24(3):955-978, 2003. 

[3] J. P. Fouque, Josselin Gamier, and A. Nachbin. Shock structure due to stochastic forcing and 
the time reversal of nonlinear waves. Physica D: Nonlinear Phenomena, 195(3-4) :324-346, 
2004. 

[4] David I. Ketcheson. High Order Strong Stability Preserving Time Integrators and Numerical 
Wave Propagation Methods for Hyperbolic PDEs. PhD thesis, Citeseer, 2009. 

[5] David I. Ketcheson and Randall J. LeVeque. Shock dynamics in layered periodic media. 
Communications in Mathematical Sciences, 10(3):859-874, 2012. 

[6] David I. Ketcheson, Kyle T. Mandli, Aron Ahmadia, Amal Alghamdi, Manuel Quezada de 
Luna, Matteo Parsani, Matthew G. Knepley, and Matthew Emmett. PyClaw: Accessible, 
Extensible, Scalable Tools for Wave Propagation Problems. SI AM Journal on Scientific 
Computing, 34(4):C210-C231, 2012. 



16 



[7] David I. Ketcheson, Matteo Parsani, and Randall J. LeVeque. High-order wave propagation 
algorithms for hyperbolic systems. SIAM Journal on Scientific Computing. In press. 

[8] Randall J. LeVeque and Marsha J Berger. Clawpack software version 4.5. 2011. Url: 
www . clawpack . org . 

[9] R.J. LeVeque. Finite volume methods for hyperbolic problems. Cambridge Univ Press, 2002. 

[10] R.J. LeVeque. Finite- volume methods for non-linear elasticity in heterogeneous media. In- 
ternational Journal for Numerical Methods in Fluids, 40(l-2):93-104, 2002. 

[11] R.J. Leveque and D.H. Yong. Solitary waves in layered nonlinear media. SIAM Journal on 
Applied Mathematics, 63(5):1539-1560, 2003. 

[12] Manuel Quezada de Luna. Nonlinear wave propagation and solitary wave formation in two- 
dimensional heterogeneous media. Master's thesis, King Abdullah University of Science and 
Technology (KAUST), 2011. 

[13] F. Santosa and W.W. Symes. A dispersive effective medium for wave propagation in periodic 
composites. SIAM Journal on Applied Mathematics, 51(4):984-1005, 1991. 

[14] N.J. Zabusky and M.D. Kruskal. Interaction of "solitons" in a collisionless plasma and the 
recurrence of initial states. Physical Review Letters, 15(6):240-243, 1965. 



17 



